!!!!#define HOMO
#define precondition
MODULE  module_gcr
!  GCR: Generalized Conjugate Residual algorithm for linear equation
!  SIAM J. NUMER. ANAL. V20,1983,pp345-357
!  equation:   Ax=f  A is nonsymmetric but with positive-definite symmetric part.
!  Method:
!  1) choose x0
!  2) compute r0=f-Ax0
!  3) Set p0=r0
!  for i=0 Step 1 Until convegence Do
!  4) ai=(ri,Api)/(Api,Api)
!  5) xi+1=xi+aipi
!  6) ri+1=ri-aiApi
!  7) compute pi+1,
!                       i  (i)
!             pi+1=ri+1+E (bj * pj)
!                      j=0
!                                        (i)
!             where  bj=-(Ari+1,Apj)/(Apj,Apj)    j=<i
!
  USE MPI
  USE h1
  USE module_halo
#define  GPTL_TIMING
#ifdef GPTL_TIMING
!   USE GPTL
#endif
  IMPLICIT NONE
include 'gptl.inc'
  !      integer, parameter,public ::r8 = selected_real_kind(13),r4 =selected_real_kind(6)
  REAL,   DIMENSION(:,:,:,:),allocatable ::   p,ap
  !REAL,   DIMENSION(:,:,:,:),ALLOCATABLE ::   a_helm
  !REAL,   DIMENSION(:,:,:),  ALLOCATABLE ::   b_helm
  !REAL,   DIMENSION(:,:,:,:),ALLOCATABLE ::   cm_helm
  !REAL,   DIMENSION(1:19,its:ite,kts-1:kte+1,jts:jte) ::   a_helm
  !REAL,   DIMENSION(its:ite,kts-1:kte+1,jts:jte)  ::   b_helm
  !REAL,   DIMENSION(its:ite,kts-1:kte+1,jts:jte,0:iter_max)  ::   ap
  !REAL,   DIMENSION(its-1:ite+1,kts-1:kte+1,jts-1:jte+1,0:iter_max)  ::   p
  !REAL,   DIMENSION(ims:ime,kms:kme,jms:jme,1:7) ::   cm_helm

CONTAINS

SUBROUTINE solve_helmholts(a_helm,b_helm,cm_helm,threshold, pi, &
                       idep, jdep,         &
                       ids,ide,jds,jde,kds,kde,   &
                       ims,ime,jms,jme,kms,kme,   &
                       its,ite,jts,jte,kts,kte,myprcid)

  IMPLICIT NONE

  INTEGER   ::           ims,ime,jms,jme,kms,kme,   &
                         ids,ide,jds,jde,kds,kde,   &
                         its,ite,jts,jte,kts,kte,   &
                         idep, jdep

  REAL*8,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(inout) :: pi
  REAL*8,DIMENSION(ims:ime,kms:kme,jms:jme):: x
!local
  REAL*4                               :: ep, threshold
  REAL*8                               :: d
  INTEGER,parameter                              :: iter_max=10
  INTEGER                              :: int1,int2
  INTEGER                              :: i,j,k,n,l,kk, myprcid
  INTEGER                              ::  ibegin,iend,jbegin,jend
  REAL,   DIMENSION(1:19,its:ite,kts-1:kte+1,jts:jte) ::   a_helm
  REAL,   DIMENSION(its:ite,kts-1:kte+1,jts:jte)  ::   b_helm
  !REAL,   DIMENSION(its:ite,kts-1:kte+1,jts:jte,0:iter_max)  ::   ap
  !REAL,   DIMENSION(its-1:ite+1,kts-1:kte+1,jts-1:jte+1,0:iter_max)  ::   p
  REAL,   DIMENSION(ims:ime,kms:kme,jms:jme,1:7) ::   cm_helm
#ifdef GPTL_TIMING
   integer   :: ret
    ret=GPTLsetoption (GPTLoverhead, 0); 
   ! ret=GPTLsetoption (GPTLdepthlimit , 2);
   !ret = gptlsetoption(gptlabort_on_error, 1)
   !ret = gptlsetutr(GPTLmpiwtime)
   ret = gptlinitialize()
   ret = gptlstart('gcr_main')
#endif

  ibegin = its
  iend   = ite
  jbegin = jts
  jend   = min (jde-1,jte)


  d=1.0d0
  ep=threshold
!    call system_clock(int1)
!    write(*,*) 'int1: ', int1
DO i=1,1
  CALL psolve_gcr_main(ep,a_helm,b_helm,cm_helm,iter_max,pi,d, &
                       idep, jdep,         &
                       ids,ide,jds,jde,kds,kde,   &
                       ims,ime,jms,jme,kms,kme,   &
                       its,ite,jts,jte,kts,kte,myprcid)
enddo
#ifdef GPTL_TIMING
   ret = gptlstop('gcr_main')
   ret = gptlpr_file("outfile_f90")
   !ret = GPTLpr_file ("outfile_f90");
   !ret = gptlpr(myprcid)
#endif


#ifdef GPTL_TIMING
   ret = gptlfinalize()
 !  WRITE(0,*) 'GRAPES-MODEL finished.'
#endif

 !    call system_clock(int2)
!    write(*,*) '(int2-int1)/1000. is ' , (int2-int1)/1000.
END SUBROUTINE solve_helmholts


SUBROUTINE initial_temp_array(iter_max,           &
                       ids,ide,jds,jde,kds,kde,   &
                       ims,ime,jms,jme,kms,kme,   &
                       its,ite,jts,jte,kts,kte)
  IMPLICIT NONE
  INTEGER,INTENT(in) :: iter_max
  INTEGER   ::           ims,ime,jms,jme,kms,kme,   &
                         ids,ide,jds,jde,kds,kde,   &
                         its,ite,jts,jte,kts,kte



  if ( .not. ALLOCATED( p ) )  ALLOCATE (p(its-1:ite+1,kts-1:kte+1,jts-1:jte+1,0:iter_max))

  if ( .not. ALLOCATED( ap ) ) ALLOCATE (ap(its:ite,kts-1:kte+1,jts:jte,0:iter_max))

  !$OMP PARALLEL WORKSHARE
  p=0.
  ap=0.
  !$OMP END PARALLEL WORKSHARE
  

END SUBROUTINE  initial_temp_array


SUBROUTINE matrixpro(a,b,c,                       &
                       idep,jdep,           &
                       ids,ide,jds,jde,kds,kde,   &
                       ims,ime,jms,jme,kms,kme,   &
                       its,ite,jts,jte,kts,kte)
  IMPLICIT NONE
  INTEGER   ::           ims,ime,jms,jme,kms,kme,   &
                         ids,ide,jds,jde,kds,kde,   &
                         its,ite,jts,jte,kts,kte, idep, jdep

  REAL,DIMENSION(1:19,its:ite,kts-1:kte+1,jts:jte),INTENT(in)  :: a
  REAL,DIMENSION(     its-1:ite+1,kts-1:kte+1,jts-1:jte+1),INTENT(inout)    :: b
  REAL,DIMENSION(     its:ite,kts-1:kte+1,jts:jte),INTENT(out)   :: c
!local
  INTEGER                     ::  ibegin,iend,jbegin,jend,i,k,j

  ibegin = its
  iend   = ite
  jbegin = jts
  jend   = min(jte,jde-1)

   !$OMP PARALLEL
   !$OMP DO PRIVATE(j,k,i) SCHEDULE(runtime) COLLAPSE(2)
  DO j=jbegin,jend
    DO k=kts,kte
      DO i=ibegin,iend
       c(i,k,j)  =                                  &
                   +a(1 ,i,k,j)*b(i  ,k  ,j  )      &
                   +a(2 ,i,k,j)*b(i-1,k  ,j  )      &
                   +a(3 ,i,k,j)*b(i+1,k  ,j  )      &
                   +a(4 ,i,k,j)*b(i  ,k  ,j-1)      &
                   +a(5 ,i,k,j)*b(i  ,k  ,j+1)      &
                   +a(6 ,i,k,j)*b(i+1,k  ,j+1)      &
                   +a(7 ,i,k,j)*b(i+1,k  ,j-1)      &
                   +a(8 ,i,k,j)*b(i-1,k  ,j-1)      &
                   +a(9 ,i,k,j)*b(i-1,k  ,j+1)      &
                   +a(10,i,k,j)*b(i  ,k-1,j  )      &
                   +a(11,i,k,j)*b(i-1,k-1,j  )      &
                   +a(12,i,k,j)*b(i+1,k-1,j  )      &
                   +a(13,i,k,j)*b(i  ,k-1,j-1)      &
                   +a(14,i,k,j)*b(i  ,k-1,j+1)      &
                   +a(15,i,k,j)*b(i  ,k+1,j  )      &
                   +a(16,i,k,j)*b(i-1,k+1,j  )      &
                   +a(17,i,k,j)*b(i+1,k+1,j  )      &
                   +a(18,i,k,j)*b(i  ,k+1,j-1)      &
                   +a(19,i,k,j)*b(i  ,k+1,j+1)
      ENDDO
    ENDDO
  ENDDO
   !$OMP END DO NOWAIT
  
    k=kts-1
   !$OMP DO PRIVATE(j,i) SCHEDULE(runtime) 
  DO j=jbegin,jend
      DO i=ibegin,iend
       c(i,k,j)  =                                  &
                   +a(1 ,i,k,j)*b(i  ,k  ,j  )      &
                   +a(2 ,i,k,j)*b(i-1,k  ,j  )      &
                   +a(3 ,i,k,j)*b(i+1,k  ,j  )      &
                   +a(4 ,i,k,j)*b(i  ,k  ,j-1)      &
                   +a(5 ,i,k,j)*b(i  ,k  ,j+1)      &
                   +a(6 ,i,k,j)*b(i+1,k  ,j+1)      &
                   +a(7 ,i,k,j)*b(i+1,k  ,j-1)      &
                   +a(8 ,i,k,j)*b(i-1,k  ,j-1)      &
                   +a(9 ,i,k,j)*b(i-1,k  ,j+1)      &
                   +a(15,i,k,j)*b(i  ,k+1,j  )      &
                   +a(16,i,k,j)*b(i-1,k+1,j  )      &
                   +a(17,i,k,j)*b(i+1,k+1,j  )      &
                   +a(18,i,k,j)*b(i  ,k+1,j-1)      &
                   +a(19,i,k,j)*b(i  ,k+1,j+1)
      ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  
    k=kte+1
  !$OMP DO PRIVATE(j,i) SCHEDULE(runtime)
  DO j=jbegin,jend
      DO i=ibegin,iend
       c(i,k,j)  =                                  &
                   +a(1 ,i,k,j)*b(i  ,k  ,j  )      &
                   +a(2 ,i,k,j)*b(i-1,k  ,j  )      &
                   +a(3 ,i,k,j)*b(i+1,k  ,j  )      &
                   +a(4 ,i,k,j)*b(i  ,k  ,j-1)      &
                   +a(5 ,i,k,j)*b(i  ,k  ,j+1)      &
                   +a(6 ,i,k,j)*b(i+1,k  ,j+1)      &
                   +a(7 ,i,k,j)*b(i+1,k  ,j-1)      &
                   +a(8 ,i,k,j)*b(i-1,k  ,j-1)      &
                   +a(9 ,i,k,j)*b(i-1,k  ,j+1)      &
                   +a(10,i,k,j)*b(i  ,k-1,j  )      &
                   +a(11,i,k,j)*b(i-1,k-1,j  )      &
                   +a(12,i,k,j)*b(i+1,k-1,j  )      &
                   +a(13,i,k,j)*b(i  ,k-1,j-1)      &
                   +a(14,i,k,j)*b(i  ,k-1,j+1)
      ENDDO
  ENDDO
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
  

END SUBROUTINE matrixpro

SUBROUTINE psolve_gcr_main(ep,a0,f0,cm_helm,iter_max,x0,d,&
                       idep, jdep,         &
                       ids,ide,jds,jde,kds,kde,   &
                       ims,ime,jms,jme,kms,kme,   &
                       its,ite,jts,jte,kts,kte,myprcid)
  IMPLICIT NONE
  INTEGER   ::           ims,ime,jms,jme,kms,kme,   &
                         ids,ide,jds,jde,kds,kde,   &
                         its,ite,jts,jte,kts,kte,   &
                         idep, jdep

  REAL*4   ,INTENT(in) :: ep
  REAL,DIMENSION(19,its:ite,kts-1:kte+1,jts:jte),INTENT(in)    :: a0
  REAL,DIMENSION(   its:ite,kts-1:kte+1,jts:jte),INTENT(in)    :: f0
  REAL*8,DIMENSION(   ims:ime,kts-1:kte+1,jms:jme),INTENT(inout) :: x0
  REAL,   DIMENSION(ims:ime,kms:kme,jms:jme,1:7) ::   cm_helm
  INTEGER,INTENT(in) :: iter_max
!  REAL,   DIMENSION(its-1:ite+1,kts-1:kte+1,jts-1:jte+1,0:iter_max)  ::   p
!  REAL,   DIMENSION(its:ite,kts-1:kte+1,jts:jte,0:iter_max)  ::   ap
! local
  REAL,DIMENSION(its:ite,kts-1:kte+1,jts:jte)               ::   r
  REAL,DIMENSION(its:ite,kts-1:kte+1,jts:jte)               ::  ar
  REAL*8,DIMENSION(0:iter_max)            ::  ac,aps,b
  REAL*8,DIMENSION(iter_max+4)            :: c1,c2
  REAL*8                                  :: d
  INTEGER                               :: ibeginh, iendh, jbeginh, jendh
  INTEGER                               :: ibegin, iend, jbegin, jend
  INTEGER                               :: i,j,k,m,l,ierr,iwb,kwb, m_tmp
  INTEGER                               :: mm, max_iteration, myprcid
  INTEGER                               :: NX,NY,NZ,NG,ILUK,NAWPP,NFWPP
  REAL*8   :: c11,c12,c1_tmp,bl,acm

  INTEGER :: jv,ret

   max_iteration = 120

   !mgrid%ims = its-1
   !mgrid%ime = ite+1
   !mgrid%jms = jts-1
   !mgrid%jme = jte+1
   !iwb=500
   !kwb=500

   ibegin = its
   iend   = ite
   jbegin = jts
   jend   = min(jte,jde-1)

   ibeginh = its-1
   iendh   = ite+1
   jbeginh = jts-1
   jendh   = min(jte+1,jde)

   NX = ite - its + 3
   NY = kte - kts + 3
   NZ = jte - jts + 3
   NG = NX * NY * NZ
   ILUK  = 0
   NAWPP = 7 * NG
   NFWPP = 23 * NG

   if ( its .ne. ids ) ibeginh = its-1
   if ( ite .ne. ide ) iendh   = ite+1
   if ( jts .ne. jds ) jbeginh = jts-1
   if ( jte .ne. jde ) jendh   = jte+1

#ifdef GPTL_TIMING
   ret=gptlstart('initialization')
#endif
  CALL initial_temp_array(iter_max,                          &
                       ids,ide,jds,jde,kds,kde,              &
                       ims,ime,jms,jme,kms,kme,              &
                       its,ite,jts,jte,kts,kte)
  p(its-1:ite+1,kts-1:kte+1,jts-1:jte+1,0) = x0(its-1:ite+1,kts-1:kte+1,jts-1:jte+1)
  m = 0

   m_tmp=m
  CALL matrixpro(a0,p(:,:,:,0),r,                                    &
                       idep, jdep,                    &
                       ids,ide,jds,jde,kds,kde,              &
                       ims,ime,jms,jme,kms,kme,              &
                       its,ite,jts,jte,kts,kte)

  d = 0.
  !$OMP PARALLEL
  !$OMP DO PRIVATE(j,k,i) SCHEDULE(runtime)
  DO j=jbegin,jend
    DO k=kts-1,kte+1
       DO i=ibegin,iend
         r(i,k,j)  = f0(i,k,j)-r(i,k,j)
       ENDDO
     ENDDO
   ENDDO
  !$OMP END DO

    !$OMP DO PRIVATE(j,k,i) SCHEDULE(runtime)
      DO j=jts,jte
        DO k=kts-1,kte+1
          DO i=its,ite
          p(i,k,j,m) = r(i,k,j)
        enddo
      enddo
    enddo
    !$OMP END DO NOWAIT
  !$OMP END PARALLEL
   
   call glob_updatehalo_real_3d(p,m,iter_max)
   m_tmp=m
   #ifdef precondition
   CALL SVRASR ( p(:,:,:,m), cm_helm,  NG, NX, NY, NZ)
   call glob_updatehalo_real_3d(p,m,iter_max)
   #endif


#ifdef GPTL_TIMING
   ret=gptlstop('initialization')
#endif

!   d = 0.
!   !$OMP PARALLEL DO PRIVATE(j,k,i) SCHEDULE(runtime)
!   DO j=jbegin,jend
!      DO k=kts-1,kte+1
!        DO i=ibegin,iend
!          !$OMP ATOMIC
!          d = d + r(i,k,j)*r(i,k,j)
!        ENDDO
!      ENDDO
!   ENDDO
!   !include "printr.inc"
!   !$OMP END PARALLEL DO
!   i=1
!   c1(1) = d
!   call MPI_ALLREDUCE(c1,c2,i,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
!   d = c2(1)
!   if(myprcid.eq.0) write(*,*)'begin of gcr', real(d,r4)

   DO mm=0,max_iteration-1
#ifdef GPTL_TIMING
   ret=gptlstart('gcr_iteration')
#endif
    m = mod(mm,iter_max-1)
    !m = mod(mm,iter_max)
    if ((m .eq. 0) .and.( mm.ne.0) ) THEN
      !$OMP PARALLEL DO PRIVATE(j,k,i) SCHEDULE(runtime)
      DO j=jts-1,jte+1
        DO k=kts-1,kte+1
          DO i=its-1,ite+1
            p(i,k,j,0) = p(i,k,j,iter_max-1)
            !p(i,k,j,0) = p(i,k,j,iter_max)
          END DO
        END DO
      END DO
      !$OMP END PARALLEL DO 
    END IF

#ifdef GPTL_TIMING
   ret=gptlstart('matrixpro')
#endif
    CALL matrixpro(a0,p(:,:,:,m),ap(:,:,:,m), &
                       idep, jdep,                    &
                       ids,ide,jds,jde,kds,kde,              &
                       ims,ime,jms,jme,kms,kme,              &
                       its,ite,jts,jte,kts,kte)

#ifdef GPTL_TIMING
   ret=gptlstop('matrixpro')
   ret=gptlstop('reduction')
#endif
    c1 = 0.
    c2 = 0.
    c11=0.d0
    c12=0.d0
    !$OMP PARALLEL DO PRIVATE(j,k,i) REDUCTION(+:c11,c12) SCHEDULE(runtime)
    DO j=jbegin,jend
      DO k=kts-1,kte+1
         DO i=ibegin,iend
           c11 = c11 + r(i,k,j)*ap(i,k,j,m)
           c12 = c12 + ap(i,k,j,m)*ap(i,k,j,m)
         END DO
      END DO
    END DO
    !$OMP END PARALLEL DO 
    
#ifdef GPTL_TIMING
   ret=gptlstop('reduction')
#endif
      c1(1)=c11
      c1(2)=c12
    i=2
    call MPI_ALLREDUCE(c1,c2,i,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)
    c1(1)=c2(1)
    c1(2)=c2(2)
    acm=c1(1)/c1(2)
    aps(m)=c1(2)

   ret=gptlstart('gcr1')
    d = 0.d0
    !$OMP PARALLEL DO PRIVATE(j,k,i) SCHEDULE(runtime)
    DO j=jbegin,jend
      DO k=kts-1,kte+1
        DO i=ibegin,iend
          x0(i,k,j) = x0(i,k,j)+acm*p(i,k,j,m)
          r(i,k,j) = r(i,k,j)-acm*ap(i,k,j,m)
        ENDDO
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO 
    !$OMP PARALLEL DO PRIVATE(j,k,i) REDUCTION(+:d) SCHEDULE(runtime)
    DO j=jbegin,jend
      DO k=kts-1,kte+1
        DO i=ibegin,iend
          d = d  + r(i,k,j)*r(i,k,j)
        ENDDO
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO 



    !$OMP PARALLEL DO PRIVATE(j,k,i) SCHEDULE(runtime)
      DO j=jts,jte
        DO k=kts-1,kte+1
          DO i=its,ite
          p(i,k,j,m+1) = r(i,k,j)
        enddo
      enddo
    enddo
    !$OMP END PARALLEL DO 
   ret=gptlstop('gcr1')
    !call glob_updatehalo(p(:,:,:,m+1),mgrid,halo=1)
   ret=gptlstart('glob_updatehalo')
   call glob_updatehalo_real_3d(p,m+1,iter_max)
   ret=gptlstop('glob_updatehalo')
   m_tmp=m+1
   !include "printp.inc"

#ifdef GPTL_TIMING
   ret=gptlstart('precondition')
#endif
   #ifdef precondition
    CALL SVRASR ( p(:,:,:,m+1), cm_helm, NG, NX, NY, NZ)
   ret=gptlstop('precondition')
   ret=gptlstart('glob_updatehalo')
   call glob_updatehalo_real_3d(p,m+1,iter_max)
   #endif

#ifdef GPTL_TIMING
   ret=gptlstop('glob_updatehalo')
#endif
   !include "printp.inc"
    !call glob_updatehalo(p(:,:,:,m+1),mgrid,halo=1)

#ifdef GPTL_TIMING
   ret=gptlstart('matrixpro')
#endif
    CALL matrixpro(a0,p(:,:,:,m+1),ar, &
                       idep, jdep,                    &
                       ids,ide,jds,jde,kds,kde,              &
                       ims,ime,jms,jme,kms,kme,              &
                       its,ite,jts,jte,kts,kte)

#ifdef GPTL_TIMING
   ret=gptlstop('matrixpro')
   ret=gptlstart('reduction')
#endif
#define MYOPT
#ifdef MYOPT
    c1 = 0.
    DO l=0,m
    c1_tmp=0.
    !$OMP PARALLEL DO PRIVATE(j,k,i) REDUCTION(+:c1_tmp) SCHEDULE(runtime)
    DO j=jbegin,jend
      DO k=kts-1,kte+1
        DO i=ibegin,iend
          c1_tmp = c1_tmp + ar(i,k,j)*ap(i,k,j,l)
        ENDDO
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO 
    c1(l+1)=c1_tmp
    ENDDO
#else
  c1 = 0.
    !$OMP PARALLEL DO PRIVATE(l,j,k,i) SCHEDULE(runtime)
    DO l=0,m
    DO j=jbegin,jend
      DO k=kts-1,kte+1
        DO i=ibegin,iend
          c1(l+1) = c1(l+1) + ar(i,k,j)*ap(i,k,j,l)
        ENDDO
      ENDDO
    ENDDO
    ENDDO
    !$OMP END PARALLEL DO 
#endif   

#ifdef GPTL_TIMING
   ret=gptlstop('reduction')
#endif
    c1(m+2) = d
    i=m+2
    call MPI_ALLREDUCE(c1,c2,i,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)
    c1=c2
    d = c2(m+2)
    !IF (  d <= ep .or.  mm == max_iteration-1 )  GOTO 10

#ifdef GPTL_TIMING
   ret=gptlstart('reductionp')
#endif
    !$OMP PARALLEL DO PRIVATE(j,l,b,k,i) SCHEDULE(runtime)
     DO j=jts-1,jte+1
      DO l=0,m
        b(l)=-c1(l+1)/aps(l)
        DO k=kts-1,kte+1
          DO i=its-1,ite+1
            p(i,k,j,m+1)=p(i,k,j,m+1)+b(l)*p(i,k,j,l)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    !$OMP END PARALLEL DO 

#ifdef GPTL_TIMING
   ret=gptlstop('reductionp')
   ret=gptlstop('gcr_iteration')
   #endif
    IF (  d <= ep .or.  mm == max_iteration-1 )  GOTO 10
  if(myprcid.eq.0)   write(*,*)'RES of gcr', real(d,r4) , 'in ',mm, 'iterations '
  ENDDO


  10   if(myprcid.eq.0)   write(*,*)'RES of gcr', real(d,r4) , 'in ',mm, 'iterations '


END SUBROUTINE psolve_gcr_main
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCC
!XXXX                                                                        XXXXC
!XXXX   AUTHOR : J. W.  CAO                                                  XXXXC
!XXXX   DATE   : Sep.11 2005                                                 XXXXC
!XXXX                                        @ Copyright BY Jianwen CAO      XXXXC
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXCC
!
!#define OPT
#ifdef OPT
      SUBROUTINE  SVRASR ( yl, b,  NG, NI, NK, NJ)
!

      IMPLICIT NONE
      INTEGER    NI,NK,NJ,NG,ILUK,NFWPP
      integer :: i,j,k,halo,l
      real :: yl(ni,nk,nj),b(ni,nk,nj,7)
   
        !$OMP PARALLEL
        !$OMP DO PRIVATE(j,i,k)
        !ni=ni-1
        !nj=nj-1
        do j=1,nj
          do i=2,ni
            yl(i,1,j)=yl(i,1,j)-b(i,1,j,2)*yl(i-1,1,j)
          enddo
           !k>1
          do k=2,nk
            do i=1,ni
              yl(i,k,j)=yl(i,k,j)-b(i,k,j,4)*yl(i,k-1,j)
            enddo

            do i=2,ni
              yl(i,k,j)=yl(i,k,j)-b(i,k,j,2)*yl(i-1,k,j)
            enddo
          enddo
        enddo
         !$OMP END DO NOWAIT
         !$OMP END PARALLEL
        !!$OMP END DO

!
          !j=1
           !k=1
           !!$OMP MASTER
           !ni=ni-l
           !nj=nj-l
           yl(ni,nk,nj)=yl(ni,nk,nj)*b(ni,nk,nj,1)
           !!$OMP END MASTER
           !!$OMP BARRIER
        !$OMP PARALLEL
           !$OMP DO PRIVATE(j,i,k)
         do j=nj,1,-1
           do i=ni-1,1,-1
             yl(i,nk,j)=(yl(i,nk,j)-b(i,nk,j,3)*yl(i+1,nk,j))*b(i,nk,j,1)
           enddo
           !k>1
           do k=nk-1,1,-1
             do i=1,ni
               yl(i,k,j)=yl(i,k,j)-b(i,k,j,5)*yl(i,k+1,j)
             enddo
             yl(ni,k,j)=yl(ni,k,j)*b(ni,k,j,1)
             do i=ni-1,1,-1
                yl(i,k,j)=(yl(i,k,j)-b(i,k,j,3)*yl(i+1,k,j))*b(i,k,j,1)
             enddo
           enddo
         enddo
         !$OMP END DO NOWAIT
         !$OMP END PARALLEL
 END SUBROUTINE  SVRASR
 #else

      SUBROUTINE  SVRASR ( yl, b,  NG, NI, NK, NJ)
!

      IMPLICIT NONE
      INTEGER    NI,NK,NJ,NG,ILUK,NFWPP
      integer :: i,j,k,halo,l
      real :: yl(ni,nk,nj),b(ni,nk,nj,7)
   
        !$OMP PARALLEL
        !$OMP DO PRIVATE(j,i,k)
        !ni=ni-1
        !nj=nj-1
        do j=1,nj
          do i=2,ni
            yl(i,1,j)=yl(i,1,j)-b(i,1,j,2)*yl(i-1,1,j)
          enddo
           !k>1
          do k=2,nk
            do i=1,ni
              yl(i,k,j)=yl(i,k,j)-b(i,k,j,4)*yl(i,k-1,j)
            enddo

            do i=2,ni
              yl(i,k,j)=yl(i,k,j)-b(i,k,j,2)*yl(i-1,k,j)
            enddo
          enddo
        enddo
        !$OMP END DO

!
          !j=1
           !k=1
           !$OMP MASTER
           !ni=ni-l
           !nj=nj-l
           yl(ni,nk,nj)=yl(ni,nk,nj)*b(ni,nk,nj,1)
           !$OMP END MASTER
           !$OMP BARRIER
           !$OMP DO PRIVATE(j,i,k)
         do j=nj,1,-1
           do i=ni-1,1,-1
             yl(i,nk,j)=(yl(i,nk,j)-b(i,nk,j,3)*yl(i+1,nk,j))*b(i,nk,j,1)
           enddo
           !k>1
           do k=nk-1,1,-1
             do i=1,ni
               yl(i,k,j)=yl(i,k,j)-b(i,k,j,5)*yl(i,k+1,j)
             enddo
             yl(ni,k,j)=yl(ni,k,j)*b(ni,k,j,1)
             do i=ni-1,1,-1
                yl(i,k,j)=(yl(i,k,j)-b(i,k,j,3)*yl(i+1,k,j))*b(i,k,j,1)
             enddo
           enddo
         enddo
         !$OMP END DO NOWAIT
         !$OMP END PARALLEL

 END SUBROUTINE  SVRASR
 #endif
! for digtial filter
END MODULE  module_GCR
